#include "base_cpp/output.h"

#include "reaction/reaction_product_enumerator.h"
#include "reaction/reaction_enumerator_state.h"

#include "base_cpp/gray_codes.h"
#include "base_c/bitarray.h"
#include "base_cpp/tlscont.h"
#include "graph/graph.h"
#include "graph/embedding_enumerator.h"
#include "graph/spanning_tree.h"
#include "molecule/molecule.h"
#include "molecule/molecule_arom.h"
#include "molecule/molecule_decomposer.h"
#include "molecule/canonical_smiles_saver.h"
#include "graph/dfs_walk.h"
#include "layout/molecule_layout.h"

ReactionEnumeratorState::ReactionMonomers::ReactionMonomers() : TL_CP_GET(_monomers), 
    TL_CP_GET(_reactant_indexes), TL_CP_GET(_deep_levels), TL_CP_GET(_tube_indexes)
{
   _monomers.clear();
   _reactant_indexes.clear();
   _deep_levels.clear();
   _tube_indexes.clear();
}

int ReactionEnumeratorState::ReactionMonomers::size()
{
   return _monomers.size();
}

const Molecule & ReactionEnumeratorState::ReactionMonomers::getMonomer( int reactant_idx, int index )
{
   int cur_idx = 0;

   for (int i = 0; i < _reactant_indexes.size(); i++)
      if (_reactant_indexes[i] == reactant_idx)
         if (cur_idx++ == index)
            return _monomers[i];

   throw Error("can't find reactant's #%d monomer #%d", reactant_idx, index);
}

const Molecule & ReactionEnumeratorState::ReactionMonomers::getMonomer( int mon_index )
{
   if (mon_index >= _monomers.size() || mon_index < 0)
      throw Error("can't find monomer #%d", mon_index);
   else
      return _monomers[mon_index];
}

void ReactionEnumeratorState::ReactionMonomers::addMonomer( int reactant_idx, 
         const Molecule &monomer, int deep_level, int tube_idx )
{
   Molecule &new_monomer = _monomers.push();

   new_monomer.cloneMolecule(monomer, NULL, NULL);
   new_monomer.calcImplicitHydrogens(true);
   new_monomer.getStereocenters().markBonds();

   _reactant_indexes.push(reactant_idx);
   _deep_levels.push(deep_level);
   _tube_indexes.push(tube_idx);
}

void ReactionEnumeratorState::ReactionMonomers::removeMonomer( int idx )
{
   for (int j = idx + 1; j < _monomers.size(); j++)
   {
      _reactant_indexes[j - 1] = _reactant_indexes[j];
      _monomers[j - 1].cloneMolecule(_monomers[j], NULL, NULL);
      _deep_levels[j - 1] = _deep_levels[j];
      _tube_indexes[j - 1] = _tube_indexes[j];
   }

   _reactant_indexes.pop();
   _monomers.pop();
   _deep_levels.pop();
   _tube_indexes.pop();
}

ReactionEnumeratorState::ReactionEnumeratorState( const Reaction &cur_reaction,
    Molecule &cur_full_product, Array<int> &cur_product_aam_array, 
    RedBlackStringMap<int> &cur_smiles_array, ReactionMonomers &cur_reaction_monomers, 
    int &cur_product_count, ObjArray< Array<int> > &cur_tubes_monomers ) : _reaction(cur_reaction),
    _smiles_array(cur_smiles_array), _reaction_monomers(cur_reaction_monomers), 
    _product_count(cur_product_count), _tubes_monomers(cur_tubes_monomers),
    TL_CP_GET(_fragments_aam_array),
    TL_CP_GET(_full_product), TL_CP_GET(_product_aam_array),
    TL_CP_GET(_product_monomers), TL_CP_GET(_is_used_hydrogens),
    TL_CP_GET(_fragments), TL_CP_GET(_is_needless_atom),
    TL_CP_GET(_att_points)
{
   _reactant_idx = _reaction.reactantBegin();
   
   _fragments_aam_array.clear();
   _full_product.clear();
   _full_product.cloneMolecule(cur_full_product, NULL, NULL);
   _fragments.clear();
   _is_needless_atom.clear();
   _att_points.clear();
   _product_aam_array.clear();
   _product_aam_array.copy(cur_product_aam_array);
   _product_monomers.clear();
   _is_used_hydrogens.clear();

   is_all_emb = false;
   is_self_react = false;
   is_one_tube = false;
   _is_frag_search = false;
   
   _deep_level = 0;
   max_deep_level = 2;
   max_product_count = 1000;

   product_proc = NULL;
   userdata = NULL;
}

ReactionEnumeratorState::ReactionEnumeratorState( ReactionEnumeratorState &cur_rpe_state ) : _reaction(cur_rpe_state._reaction), 
    _smiles_array(cur_rpe_state._smiles_array), 
    _reaction_monomers(cur_rpe_state._reaction_monomers),
    _product_count(cur_rpe_state._product_count),
    _tubes_monomers(cur_rpe_state._tubes_monomers),
    TL_CP_GET(_fragments_aam_array),
    TL_CP_GET(_full_product), TL_CP_GET(_product_aam_array),
    TL_CP_GET(_product_monomers), TL_CP_GET(_is_used_hydrogens),
    TL_CP_GET(_fragments), TL_CP_GET(_is_needless_atom), 
    TL_CP_GET(_att_points)
{
   _reactant_idx = cur_rpe_state._reactant_idx;
   
   _fragments_aam_array.clear();
   _fragments_aam_array.copy(cur_rpe_state._fragments_aam_array);
   _full_product.clear();
   _full_product.cloneMolecule(cur_rpe_state._full_product, NULL, NULL);
   _fragments.clear();
   _fragments.cloneMolecule(cur_rpe_state._fragments, NULL, NULL);
   _is_needless_atom.clear();
   _is_needless_atom.copy(cur_rpe_state._is_needless_atom);
   _att_points.clear();
   for (int i = 0; i < cur_rpe_state._att_points.size(); i++)
   {
      Array<int> &new_array = _att_points.push();
      new_array.copy(cur_rpe_state._att_points[i]);
   }
   _product_aam_array.clear();
   _product_aam_array.copy(cur_rpe_state._product_aam_array);
   _product_monomers.clear();
   _product_monomers.copy(cur_rpe_state._product_monomers);
    _is_used_hydrogens.clear();
   _is_used_hydrogens.copy(cur_rpe_state._is_used_hydrogens);
   
   max_product_count = cur_rpe_state.max_product_count;
   _deep_level = cur_rpe_state._deep_level;
   max_deep_level = cur_rpe_state.max_deep_level;
   is_all_emb = cur_rpe_state.is_all_emb;
   is_self_react = cur_rpe_state.is_self_react;
   is_one_tube = cur_rpe_state.is_one_tube;

   _is_frag_search = false;

   product_proc = cur_rpe_state.product_proc;
   userdata = cur_rpe_state.userdata;
}

int ReactionEnumeratorState::buildProduct( void )
{
   int cnt = 0;

   if (_product_count >= max_product_count)
      return 0;

   if (_reactant_idx == _reaction.reactantEnd())
   {
      /* Product is ready */
      _productProcess();
      int max_depth = 0;
      return 0;
   }

   int monomer_idx = 0;
   for (int i = 0; i < _reaction_monomers._monomers.size(); i++)
   {
      if (!is_one_tube)
      {
         int j;
         int tube_idx = -1;
         for (j = 0; j < _product_monomers.size(); j++)
         {
            if (_reaction_monomers._reactant_indexes[_product_monomers[j]] == 
                _reaction_monomers._reactant_indexes[i])
                break;
         
            if (tube_idx == -1)
               tube_idx = _reaction_monomers._tube_indexes[_product_monomers[j]];
         }

         if (tube_idx != -1)
            if (_tubes_monomers[tube_idx].find(i) == -1)
               continue;
         
         if (j != _product_monomers.size())
            continue;
      }

      if (!is_self_react)
         if ((_reaction_monomers._deep_levels[i] == 0) && (_product_monomers.find(i) != -1))
            continue;

      ReactionEnumeratorState rpe_state(*this);

      rpe_state._deep_level += _reaction_monomers._deep_levels[i];

      if (rpe_state._deep_level > rpe_state.max_deep_level)
         return 0;

      rpe_state._product_monomers.push(i);

      QS_DEF(Molecule, ee_monomer);
      ee_monomer.clear();
      ee_monomer.cloneMolecule(_reaction_monomers._monomers[i], NULL, NULL);
      ee_monomer.calcImplicitHydrogens(true);
      ee_monomer.cis_trans.build(ee_monomer, NULL);
      //ee_monomer.unfoldHydrogens(NULL);

      rpe_state._start_ee(ee_monomer);
   }

   return 0;
}

void ReactionEnumeratorState::_productProcess( void )
{
   if (_deep_level + 1 >= max_deep_level)
      return;

   QS_DEF(Molecule, ready_product);
   ready_product.clear();
   if (!_attachFragments(ready_product))
      return;

   _foldHydrogens(ready_product);

   MoleculeLayout mol_layout(ready_product);
   mol_layout.respect_existing_layout = false;
   mol_layout.bond_length = 2;
   mol_layout.make();
   ready_product.getStereocenters().markBonds();

   /* Adding a product to monomers lists */
   if (is_all_emb)
   {
      int tube_idx = -1;
      for (int i = 0; i < _product_monomers.size(); i++)
      {
         int used_mon = _product_monomers[i];
         if ((tube_idx = _reaction_monomers._tube_indexes[used_mon]) != -1)
            break;
      }
      if (tube_idx == -1)
      {
         tube_idx = _tubes_monomers.size();
         _tubes_monomers.push().copy(_product_monomers);
      }

      for (int i = _reaction.reactantBegin(); i != _reaction.reactantEnd(); i = _reaction.reactantNext(i))
      {
         _tubes_monomers[tube_idx].push(_reaction_monomers.size());
         _reaction_monomers.addMonomer(i, ready_product, _deep_level + 1, tube_idx);
      }
   }

   QS_DEF(Array<char>, cur_smiles);
   cur_smiles.clear();
   {
      ArrayOutput arr_out(cur_smiles);
      CanonicalSmilesSaver product_cs_saver(arr_out);
      product_cs_saver.saveMolecule(ready_product);
   }
   cur_smiles.push(0);
   if (_smiles_array.find(cur_smiles.ptr()))
   {
      int *found_count = _smiles_array.at2(cur_smiles.ptr());
      *found_count++;
      return;
   }
   _product_count++;
   _smiles_array.insert(cur_smiles.ptr(), 1);

   if (product_proc != NULL)
      product_proc(ready_product, _product_monomers, userdata);
}

void ReactionEnumeratorState::_foldHydrogens( Molecule &molecule )
{
   for (int i = molecule.vertexBegin(); i != molecule.vertexEnd(); i = molecule.vertexNext(i))
   {
      const Atom &atom = molecule.getAtom(i);

      if (atom.pureHydrogen())
         molecule.removeAtom(i);
   }
   molecule.calcImplicitHydrogens(true);
}

bool ReactionEnumeratorState::_nextMatchProcess( EmbeddingEnumerator &ee, const Molecule &monomer )
{
   ReactionEnumeratorState rpe_state(*this);
  
   rpe_state._is_used_hydrogens.clear_resize(monomer.vertexEnd());
   rpe_state._is_used_hydrogens.zerofill();
   rpe_state._is_frag_search = _is_frag_search;

   ee.userdata = &rpe_state;

   bool stop_flag = ee.processNext();

   return stop_flag;
}

void ReactionEnumeratorState::_markAromaticBonds( Molecule &molecule )
{
   AromatizerBase arom_base(molecule);
   arom_base.aromatize();

   for (int i = molecule.edgeBegin(); i != molecule.edgeEnd(); i = molecule.edgeNext(i))
      if (arom_base.isBondAromatic(i))
         molecule.setBondAromatic(i);
}

void ReactionEnumeratorState::_start_ee( const Molecule &monomer )
{
   QS_DEF(Molecule, ee_reactant);
   ee_reactant.clear();
   ee_reactant.cloneMolecule(_reaction.getMolecule(_reactant_idx), NULL, NULL);
   ee_reactant.calcImplicitHydrogens(true);
   ee_reactant.cis_trans.build(ee_reactant, NULL);
   //ee_reactant.unfoldHydrogens(NULL);
   ee_reactant.findBondsInRings();

   /* Finding in reactant query atoms with two neighbors */
   QS_DEF(Array<int>, qa_array);
   qa_array.clear();
   for (int j = ee_reactant.vertexBegin(); j != ee_reactant.vertexEnd(); j = ee_reactant.vertexNext(j))
   {
      const QueryAtom &qa = ee_reactant.getQueryAtom(j);
      const Vertex &vertex = ee_reactant.getVertex(j);

      if (qa.type != QUERY_ATOM_RGROUP)
         continue;

      if (vertex.degree() > 2)
         throw Error("query atom can't have more than two neighbors");

      if (vertex.degree() == 2)
         qa_array.push(j);
   }

   GrayCodesEnumerator gce(qa_array.size(), true);

   /* Enumerating of different any atoms existence variants */
   while (!gce.isDone())
   {
      int changed_idx = gce.getBitChangeIndex();

      if (changed_idx != -1)
      {
         const byte *code = gce.getCode();
         int is_double = bitGetBit(code, changed_idx);
         int change_atom_idx = qa_array[changed_idx];
         
         _changeQueryNode(ee_reactant, change_atom_idx, is_double);
      }

      //_markAromaticBonds(ee_monomer);
      //_markAromaticBonds(ee_reactant);

      QS_DEF(Molecule, ee_monomer);
      ee_monomer.clear();

      ee_monomer.cloneMolecule(monomer, NULL, NULL);
      ee_monomer.calcImplicitHydrogens(true);

      ee_monomer.findBondsInRings();

      EmbeddingEnumerator ee(ee_monomer);
      ee.cb_embedding = _embeddingCallback;
      ee.cb_match_vertex = _matchVertexCallback;
      ee.cb_match_edge = _matchEdgeCallback;
      ee.setSubgraph(ee_reactant);
      
      ee.processStart();

      while (true)
      {
         bool stop_flag = _nextMatchProcess(ee, ee_monomer);
         gce.next();
         if (!stop_flag)
            break;
      }
   }
}

void ReactionEnumeratorState::_changeQueryNode( Molecule &ee_reactant, int change_atom_idx, int is_double )
{
   QS_DEF(Molecule, reactant_copy);
   reactant_copy.clear();
   reactant_copy.cloneMolecule(ee_reactant, NULL, NULL);

   if (is_double)
   {
      /* Making one query atom for one bond */
      const Vertex &vertex = ee_reactant.getVertex(change_atom_idx);
      const Atom &atom = ee_reactant.getAtom(change_atom_idx);
      
      for (int j = vertex.neiNext(vertex.neiBegin()); j != vertex.neiEnd(); j = vertex.neiNext(j))
      {
         int nv_idx = vertex.neiVertex(j);
         
         MoleculeRGroups *rgroups = reactant_copy.getRGroups();
         const Array<int> &rg_array = rgroups->getSiteRGroups(change_atom_idx);

         int new_atom_idx = reactant_copy.addAtomCopy(atom);
         QueryAtom &new_qa = reactant_copy.getQueryAtom2(new_atom_idx);
         new_qa.type = QUERY_ATOM_RGROUP;
         rgroups->addRGroupOnSite(new_atom_idx, rg_array[0]);

         int old_bond_idx = reactant_copy.findEdgeIndex(nv_idx, change_atom_idx);
         bool old_bond_in_ring = reactant_copy.bondInRing(old_bond_idx);

         reactant_copy.flipBond(nv_idx, change_atom_idx, new_atom_idx);

         int new_bond_idx = reactant_copy.findEdgeIndex(nv_idx, new_atom_idx);
         reactant_copy.setBondInRing(new_bond_idx, old_bond_in_ring);
      }
   }
   else
   {
      /* Making one query atom for all neighbor bonds */
      const Array<int> &rg_array = ee_reactant.getRGroups()->getSiteRGroups(change_atom_idx);
                    
      for (int j = ee_reactant.vertexBegin(); j != ee_reactant.vertexEnd(); 
               j = ee_reactant.vertexNext(j))
      {
         const QueryAtom &qa = ee_reactant.getQueryAtom(j);

         if (qa.type != QUERY_ATOM_RGROUP)
            continue;

         MoleculeRGroups *rgroups = reactant_copy.getRGroups();
         const Array<int> &cur_rg_array = rgroups->getSiteRGroups(j);
         
         if ((j != change_atom_idx) && (cur_rg_array[0] == rg_array[0]))
         {
            const Vertex &vertex = reactant_copy.getVertex(j);
            int nv_idx = vertex.neiVertex(vertex.neiBegin());
            int ne_idx = vertex.neiEdge(vertex.neiBegin());
            const Bond &bond = reactant_copy.getBond(ne_idx);

            reactant_copy.flipBond(nv_idx, j, change_atom_idx);
            reactant_copy.removeAtom(j);
         }
      }
   }
   ee_reactant.cloneMolecule(reactant_copy, NULL, NULL);
}

bool ReactionEnumeratorState::_matchVertexCallback( const Graph &subgraph, const Graph &supergraph,
                                                 const int *core_sub, int sub_idx, int super_idx, void *userdata )
{
   const ReactionEnumeratorState *rpe_state = (ReactionEnumeratorState *)userdata;
   const Molecule &supermolecule = (const Molecule &)supergraph;
   const Molecule &submolecule = (const Molecule &)subgraph;
   const Atom &a_super = supermolecule.getAtom(super_idx);
   const Vertex &v_super = supermolecule.getVertex(super_idx);
   const Atom &a_sub = submolecule.getAtom(sub_idx);
   const Vertex &v_sub = submolecule.getVertex(sub_idx);
   const QueryAtom &qa_sub = submolecule.getQueryAtom(sub_idx);

   if (submolecule.getStereocenters().getType(sub_idx) > supermolecule.getStereocenters().getType(super_idx))
      return false;

   if (supermolecule.getAtom(super_idx).pureHydrogen() && 
        ((submolecule.getAtom(sub_idx).pureHydrogen()) || (qa_sub.type == QUERY_ATOM_RGROUP) || 
         (qa_sub.type == QUERY_ATOM_A) || 
         ((qa_sub.type == QUERY_ATOM_LIST) && (qa_sub.listContains(a_super.label)))))
   {
      if (v_super.degree() != 1)
         throw Error("Hydrogen have more than one neighbor");
   
      const Vertex &v_super_nei = supermolecule.getVertex(v_super.neiVertex(v_super.neiBegin()));

      for (int i = v_super_nei.neiBegin(); i != v_super_nei.neiEnd(); i = v_super_nei.neiNext(i))
      {
         int nv_idx = v_super_nei.neiVertex(i);

         if (!supermolecule.getAtom(nv_idx).pureHydrogen())
            continue;

         if (rpe_state->_is_used_hydrogens[nv_idx])
            continue;
         if (nv_idx == super_idx)
         {
            rpe_state->_is_used_hydrogens[super_idx] = 1;
            break;
         }

         return false;
      }
   }

   if (qa_sub.type == QUERY_ATOM_RGROUP)
      return true;

   if (v_super.degree() != v_sub.degree())
      return false;

   /* We are checking already founded atom */
   if (rpe_state->_is_frag_search)
      if (rpe_state->_is_needless_atom[super_idx])
         return false;

   if ((qa_sub.type == QUERY_ATOM_A) || 
      ((qa_sub.type == QUERY_ATOM_LIST) && (qa_sub.listContains(a_super.label))))
      return true;

   if (a_super.label != a_sub.label)
      return false;
   if (a_super.charge != a_sub.charge)
      return false;

   return true;
}

bool ReactionEnumeratorState::_matchEdgeCallback( const Graph &subgraph, const Graph &supergraph,
                                               int self_idx, int other_idx, void *userdata )
{
   const ReactionEnumeratorState *rpe_state = (ReactionEnumeratorState *)userdata;
   const Molecule &supermolecule = (const Molecule &)supergraph;
   const Molecule &submolecule = (const Molecule &)subgraph;
   const Bond &b_other = supermolecule.getBond(other_idx);
   const Bond &b_self = submolecule.getBond(self_idx);
   const QueryBond &qb_other = supermolecule.getQueryBond(other_idx);
   const QueryBond &qb_self = submolecule.getQueryBond(self_idx);
   
   if (supermolecule.bondInRing(other_idx) != submolecule.bondInRing(self_idx))
      return false;

   if (qb_self.type == QUERY_BOND_ANY)
      return true;

   if ((b_self.type == BOND_AROMATIC) && (b_other.type == BOND_AROMATIC))
      return true;
   
   if (b_self.type != b_other.type)
      return false;

   return true;
}

void ReactionEnumeratorState::_findFragAtoms( Array<byte> &unfrag_mon_atoms, 
                   const Molecule &submolecule, int *core_sub, int *core_super )
{
   for (int i = submolecule.vertexBegin(); i != submolecule.vertexEnd(); i = submolecule.vertexNext(i))
   {
      const QueryAtom &qa = submolecule.getQueryAtom(i);

      if (qa.type != QUERY_ATOM_RGROUP)
         unfrag_mon_atoms[core_sub[i]] = 1;
   }
}

bool ReactionEnumeratorState::_isFragmentCorrect( const MoleculeRGroups *sub_rgroups, 
                   const MoleculeDecomposer &frag_decomp, const Array<int> &sub_qa_array, 
                   const int *mapping )
{
   QS_DEF(Array<int>, used_components);
   used_components.clear_resize(frag_decomp.getComponentsCount());
   used_components.fffill();

   for (int i = 0; i < sub_qa_array.size(); i++)
   {
      const Array<int> &rg_array = sub_rgroups->getSiteRGroups(sub_qa_array[i]);
      int component_idx = frag_decomp.getComponent(mapping[sub_qa_array[i]]);

      if (used_components[component_idx] == -1)
         used_components[component_idx] = rg_array[0];
      else if (used_components[component_idx] != rg_array[0])
         return false;
   }

   return true;
}

int ReactionEnumeratorState::_findAttachmentPoint( int idx )
{
   for (int i = 0; i < _att_points.size(); i++)
   {
      int rg_idx = _att_points[i].find(idx);
      if (rg_idx != -1)
         return i;
   }
   return -1;
}

int ReactionEnumeratorState::_findRGroupsReactant( int rg_idx )
{
   for (int j = _reaction.reactantBegin(); j != _reaction.reactantEnd(); j = _reaction.reactantNext(j))
   {
      const Molecule &reactant = _reaction.getMolecule(j);

      for (int k = reactant.vertexBegin(); k != reactant.vertexEnd(); k = reactant.vertexNext(k))
      {
         if (reactant.getQueryAtom(k).type != QUERY_ATOM_RGROUP)
            continue;

         if (reactant.getRGroups()->getSiteRGroups(k)[0] == rg_idx)
            return j;
      }
   }

   return -1;
}

void ReactionEnumeratorState::_cleanFragments( void )
{
   QS_DEF(Array<int>, needless_bonds);
   needless_bonds.clear();
   
   for (int i = _fragments.vertexBegin(); i != _fragments.vertexEnd(); i = _fragments.vertexNext(i))
   {
      if (_is_needless_atom[i])
      {
         _fragments.removeAtom(i);
         continue;
      }

      /* if we have R1-R2 variant bond should be removed (if R1 and R2 are in one reactant) */
      int rg_idx1 = _findAttachmentPoint(i);
      if (rg_idx1 == -1)
         continue;

      int reactant1 = _findRGroupsReactant(rg_idx1);

      const Vertex &v = _fragments.getVertex(i);
      for (int j = v.neiBegin(); j != v.neiEnd(); j = v.neiNext(j))
      {
         int nv_idx = v.neiVertex(j);
         int ne_idx = v.neiEdge(j);

         if (needless_bonds.find(ne_idx) != -1)
            continue;

         int rg_idx2 = _findAttachmentPoint(nv_idx);
         if (rg_idx2 != -1)
         {
            int reactant2 = _findRGroupsReactant(rg_idx2);

            if (reactant1 == reactant2)
               needless_bonds.push(ne_idx);
         }
      }
   }

   for (int j = 0; j < needless_bonds.size(); j++)
      _fragments.removeBond(needless_bonds[j]);
}

void ReactionEnumeratorState::_findR2PMapping( const Molecule &reactant, const  Molecule &product, Array<int> &mapping )
{
   const Array<int> &reactant_aam_array = _reaction.getAAMArray(_reactant_idx);

   for(int i = reactant.vertexBegin(); i != reactant.vertexEnd(); i = reactant.vertexNext(i))
   {
      if ((i >= reactant_aam_array.size()) && (reactant.getQueryAtom(i).type != QUERY_ATOM_RGROUP))
         break;

      if (reactant.getQueryAtom(i).type == QUERY_ATOM_RGROUP)
      {
         for (int j = product.vertexBegin(); j != product.vertexEnd(); j = product.vertexNext(j))
         {
            if (product.getQueryAtom(j).type != QUERY_ATOM_RGROUP)
               continue;

            int pr_rg_idx = product.getRGroups()->getSiteRGroups(j)[0];
            int sub_rg_idx = reactant.getRGroups()->getSiteRGroups(i)[0];
            if (pr_rg_idx != sub_rg_idx)
               continue;

            mapping[i] = j;
            break;
         }
      }
      else if (reactant_aam_array[i] != 0)
         mapping[i] = _product_aam_array.find(reactant_aam_array[i]);
   }
}

void ReactionEnumeratorState::_invertStereocenters( Molecule &molecule, int edge_idx )
{
   const Edge &edge = molecule.getEdge(edge_idx);
   int edge_end_idx = edge.end;
   const Vertex &edge_end = molecule.getVertex(edge_end_idx);
   int other_end_idx = edge.findOtherEnd(edge_end_idx);

   QS_DEF(Array<int>, was_atoms);
   was_atoms.clear_resize(molecule.vertexEnd());
   was_atoms.zerofill();

   for (int i = edge_end.neiBegin(); i != edge_end.neiEnd(); i = edge_end.neiNext(i))
   {
      int nei_atom_idx = edge_end.neiVertex(i);

      if (nei_atom_idx == other_end_idx)
         continue;

      QS_DEF(Array<int>, ignored_atoms);
      ignored_atoms.clear_resize(molecule.vertexEnd());
      ignored_atoms.zerofill();
      ignored_atoms[edge_end_idx] = 1;
      QS_DEF(Array<int>, atom_ranks);
      atom_ranks.clear_resize(molecule.vertexEnd());
      atom_ranks.zerofill();
      atom_ranks[nei_atom_idx] = -1;

      DfsWalk dfs(molecule);
      dfs.ignored_vertices = ignored_atoms.ptr();
      dfs.vertex_ranks = atom_ranks.ptr();
      dfs.walk();
      const Array<DfsWalk::SeqElem> &atoms_to_reflect = dfs.getSequence();

      for (int j = 0; j < atoms_to_reflect.size(); j++)
      {
         if (((j > 0) && (atoms_to_reflect[j].parent_vertex == -1)) || 
             was_atoms[atoms_to_reflect[j].idx])
            break;

         if (molecule.getStereocenters().exists(atoms_to_reflect[j].idx))
            molecule.getStereocenters().inversePyramid(atoms_to_reflect[j].idx);

         was_atoms[atoms_to_reflect[j].idx] = 1; 
      }
   }
}

void ReactionEnumeratorState::_cistransUpdate( const Molecule &submolecule,
                                Molecule &supermolecule, int *frag_mapping, 
                                const Array<int> &rp_mapping, int *core_sub)
{
   QS_DEF(Array<int>, cistrans_changed_bonds);
   cistrans_changed_bonds.clear();

   for(int i = submolecule.edgeBegin(); i != submolecule.edgeEnd(); i = submolecule.edgeNext(i))
   {
      if (!MoleculeCisTrans::isGeomStereoBond(submolecule, i, NULL, false))
         continue;

      const Edge &edge = submolecule.getEdge(i);
      const int *subs = submolecule.cis_trans.getSubstituents(i);

      if ((rp_mapping[edge.beg] == -1) || (rp_mapping[edge.end] == -1))
      {
         continue;
         //or throw Error("Incorrect AAM on stereo bond");
      }
      for (int j = 0; j < 2; j++)
      {
        if ((subs[j] != -1) && (rp_mapping[subs[j]] != -1) && 
            (_full_product.findEdgeIndex(rp_mapping[subs[j]], rp_mapping[edge.beg]) == -1))
           return;
        if ((subs[j + 2] != -1) && (rp_mapping[subs[j + 2]] != -1) && 
            (_full_product.findEdgeIndex(rp_mapping[subs[j + 2]], rp_mapping[edge.end]) == -1))
           return;
      }
      int ss_sign, sp_sign;

      ss_sign = MoleculeCisTrans::getMappingParitySign(submolecule, supermolecule, i, core_sub);
      sp_sign = MoleculeCisTrans::getMappingParitySign(submolecule, _full_product, i, rp_mapping.ptr());

      if (sp_sign > 0)
         continue;

      int product_edge_idx = _full_product.findMappedEdge(submolecule, _full_product, i, rp_mapping.ptr());
      if (product_edge_idx == -1)
         continue;

      if (_full_product.getQueryBond(product_edge_idx).stereo_care)
         continue;

      if (sp_sign * ss_sign > 0)
      {
         int new_parity = _full_product.cis_trans.getParity(product_edge_idx);
         new_parity = (new_parity == MoleculeCisTrans::CIS ? MoleculeCisTrans::TRANS : MoleculeCisTrans::CIS);
         _full_product.cis_trans.setParity(product_edge_idx, new_parity);

         int super_edge_idx = supermolecule.findMappedEdge(submolecule, supermolecule, i, core_sub);
         _invertStereocenters(supermolecule, super_edge_idx);
      }
   }
}

void ReactionEnumeratorState::_stereocentersUpdate( const Molecule &submolecule,
                   const Molecule &supermolecule, const Array<int> &rp_mapping,
                   int *core_sub, int *core_super )
{
   QS_DEF(Array<int>, mp_mapping);
   mp_mapping.clear_resize(supermolecule.vertexEnd());
   mp_mapping.fffill();

   /* Finding of monomer to product atom to atom mapping */
   for(int i = supermolecule.vertexBegin(); i != supermolecule.vertexEnd(); i = supermolecule.vertexNext(i))
      mp_mapping[i] = ((core_super[i] != -1) ? rp_mapping[core_super[i]] : -1);

   for (int i = submolecule.vertexBegin(); i != submolecule.vertexEnd(); i = submolecule.vertexNext(i))
   {
      int sub_ex = 1, pr_ex = 1;

      if (submolecule.getQueryAtom(i).type == QUERY_ATOM_RGROUP)
         continue;

      if (!submolecule.getStereocenters().exists(i))
         sub_ex = -1;
      if (!_full_product.getStereocenters().exists(rp_mapping[i]))
         pr_ex = -1;
      if (!supermolecule.getStereocenters().exists(core_sub[i]))
         continue;

      int mon_type, mon_group, mon_pyramid[4];
      supermolecule.getStereocenters().get(core_sub[i], mon_type, mon_group, mon_pyramid);

      int new_pr_pyramid[4];
      for (int j = 0; j < 4; j++)
         new_pr_pyramid[j] = ((mon_pyramid[j] != -1) ? mp_mapping[mon_pyramid[j]] : -1);

      if (sub_ex * pr_ex < 0)
      {
         continue;
      }
      else if ((sub_ex < 0) && (pr_ex < 0))
      {
         /* if there is no stereo info in reaction take monomer stereo info*/
         _full_product.getStereocenters().add(rp_mapping[i], mon_type, mon_group, new_pr_pyramid);
         continue;
      }
      else
      {
         int rct_type = submolecule.getStereocenters().getType(i);
         int pr_type = _full_product.getStereocenters().getType(rp_mapping[i]);
         int pr_group = _full_product.getStereocenters().getGroup(rp_mapping[i]);

         if ((rct_type == MoleculeStereocenters::ATOM_ANY) || (pr_type == MoleculeStereocenters::ATOM_ANY))
            continue;
         if (pr_type == MoleculeStereocenters::ATOM_ABS)
            continue;
         if (pr_type != MoleculeStereocenters::ATOM_OR)
         {
            pr_type = mon_type;
            pr_group = mon_group;
         }

         int mapping[4];
         /* Reactant to product stereocenter's pyramid mapping finding */
         MoleculeStereocenters::getPyramidMapping(submolecule.getStereocenters(),
            _full_product.getStereocenters(), i, rp_mapping.ptr(), mapping, false);

         _full_product.getStereocenters().remove(rp_mapping[i]);
         _full_product.getStereocenters().add(rp_mapping[i], pr_type, pr_group, new_pr_pyramid);
         
         if (MoleculeStereocenters::isPyramidMappingRigid(mapping))
            continue;

         _full_product.getStereocenters().inversePyramid(rp_mapping[i]);
      }
   }
}

void ReactionEnumeratorState::_anyReplace( Molecule &product, Molecule &uncleaned_fragments )
{
   /* ANY and LIST query atoms replacing */
   for (int i = product.vertexBegin(); i != product.vertexEnd(); i = product.vertexNext(i))
   {
      const QueryAtom &pr_qa = product.getQueryAtom(i);
      Atom &pr_a = product.getAtom2(i);

      if ((pr_qa.type != QUERY_ATOM_A) && (pr_qa.type != QUERY_ATOM_LIST))
         continue;

      int pr_qa_aam = _product_aam_array[i];
      if (pr_qa_aam == 0)
         throw Error("Incorrect AAM");
      int frags_qa_idx = _fragments_aam_array.find(pr_qa_aam);
      if (frags_qa_idx == -1)
         throw Error("Incorrect AAM");

      pr_a.copy(uncleaned_fragments.getAtom(frags_qa_idx));
      product.getQueryAtom2(i).type = 0;
   }

   /* ANY qwery bonds replacing */
   for (int i = product.edgeBegin(); i != product.edgeEnd(); i = product.edgeNext(i))
   {
      const QueryBond &pr_qb = product.getQueryBond(i);
      Bond &pr_b = product.getBond2(i);
      const Edge &pr_edge = product.getEdge(i);

      if (pr_qb.type != QUERY_BOND_ANY)
         continue;

      int pr_beg_aam = _product_aam_array[pr_edge.beg];
      int pr_end_aam = _product_aam_array[pr_edge.end];
      if ((pr_beg_aam == 0) || (pr_end_aam == 0))
         throw Error("Incorrect AAM");
      int frags_beg = _fragments_aam_array.find(pr_beg_aam);
      int frags_end = _fragments_aam_array.find(pr_end_aam);
      if ((frags_beg == -1) || (frags_end == -1))
         throw Error("Incorrect AAM");
      int frags_bond_idx = uncleaned_fragments.findEdgeIndex(frags_beg, frags_end);
      if (frags_bond_idx == -1)
         throw Error("Incorrect AAM");
      
      /* Here should be atom copying */
      pr_b.type = uncleaned_fragments.getBond(frags_bond_idx).type;
      product.getQueryBond2(i).type = 0;
   }
}

void ReactionEnumeratorState::_completeCisTranse( Molecule &product, Molecule &uncleaned_fragments, 
                                                  Array<int> &frags_mapping )
{
   for (int i = _fragments.edgeBegin(); i != _fragments.edgeEnd(); i = _fragments.edgeNext(i))
   {
      if ((_fragments.getBond(i).type != BOND_DOUBLE) ||
         (_fragments.cis_trans.getParity(i) != 0) ||
         (uncleaned_fragments.cis_trans.getParity(i) == 0))
         continue;

      const Edge &edge = uncleaned_fragments.getEdge(i);
      int super_parity = uncleaned_fragments.cis_trans.getParity(i);
      const int *subs = uncleaned_fragments.cis_trans.getSubstituents(i);
      int new_subs[4];

      int k;
      for (k = 0; k < 4; k++)
      {
         if (subs[k] == -1)
         {
            new_subs[k] = -1;
            continue;
         }

         if (frags_mapping[subs[k]] == -1)
         {
            int sub_aam = _fragments_aam_array[subs[k]];
            if (sub_aam == 0)
               break;
            int substiuent = _product_aam_array.find(sub_aam);
            if (substiuent == -1)
               break;
            new_subs[k] = substiuent;
            continue;
         }
         new_subs[k] = frags_mapping[subs[k]];
      }

      if (k < 4)
         continue;

      int pr_bond_idx = product.findEdgeIndex(frags_mapping[edge.beg], frags_mapping[edge.end]);
      product.cis_trans.add(pr_bond_idx, new_subs, uncleaned_fragments.cis_trans.getParity(i));
   }
}

bool ReactionEnumeratorState::_attachFragments( Molecule &product )
{
   product.clear();
   product.cloneMolecule(_full_product, NULL, NULL);

   QS_DEF(Molecule, uncleaned_fragments);
   uncleaned_fragments.clear();
   uncleaned_fragments.cloneMolecule(_fragments, NULL, NULL);

   /* ANY and LIST query atoms and bonds replacing */
   _anyReplace(product, uncleaned_fragments);

   _cleanFragments();

   QS_DEF(Array<int>, frags_mapping);
   frags_mapping.clear_resize(_fragments.vertexEnd());
   frags_mapping.fffill();
   product.mergeWithMolecule(_fragments, &frags_mapping);

   for (int i = _full_product.vertexBegin(); i != _full_product.vertexEnd(); i = _full_product.vertexNext(i))
   {
      const QueryAtom &pr_qa = _full_product.getQueryAtom(i);
      if (pr_qa.type == QUERY_ATOM_RGROUP)
      {
         int pr_rgroup = _full_product.getRGroups()->getSiteRGroups(i)[0];
         const Vertex &pr_v = product.getVertex(i);
         int nei_idx = pr_v.neiBegin();

         QS_DEF(Array<int>, pr_neibours);
         pr_neibours.clear();
         for (int j = pr_v.neiBegin(); j != pr_v.neiEnd(); j = pr_v.neiNext(j))
            pr_neibours.push(pr_v.neiVertex(j));

         if (pr_neibours.size() > 2)
            throw Error("RGroup atom can't have more than two neibours");

         /* Setting the order of rgroup atom neibours by AAM (less is first) */
         if (pr_neibours.size() == 2)
         {
            if (((_product_aam_array.size() > pr_neibours[0]) && 
                (_product_aam_array.size() > pr_neibours[1])) &&
                _product_aam_array[pr_neibours[0]] > _product_aam_array[pr_neibours[1]])
            {
               int tmp = pr_neibours[0];
               pr_neibours[0] = pr_neibours[1];
               pr_neibours[1] = tmp;
            }
         }
         
         for (int j = 0; j < pr_neibours.size(); j++)
         {
            if (product.findEdgeIndex(pr_neibours[j], frags_mapping[_att_points[pr_rgroup][j]]) != -1)
               return false;

            product.flipBond(pr_neibours[j], i, frags_mapping[_att_points[pr_rgroup][j]]);
         }
         product.removeAtom(i);

         int nv_idx = 0;
         for (int j = 0; j < _att_points[pr_rgroup].size(); j++)
         {
            if (uncleaned_fragments.getStereocenters().exists(_att_points[pr_rgroup][j]))
            {
               int type, group, pyramid[4];
               uncleaned_fragments.getStereocenters().get(_att_points[pr_rgroup][j], type, group, pyramid);

               int new_pyramid[4];

               for (int k = 0; k < 4; k++)
               {
                  if (pyramid[k] == -1)
                     new_pyramid[k] = -1;
                  else if (!_is_needless_atom[pyramid[k]])
                     new_pyramid[k] = frags_mapping[pyramid[k]];
                  else
                     new_pyramid[k] = pr_neibours[nv_idx++];
               }

               product.getStereocenters().add(frags_mapping[_att_points[pr_rgroup][j]], type, group, new_pyramid);
            }

            if (nv_idx == 2)
               break;
         }
      }
   }

   /* Updating of cis-trans information on product & monomer's fragment border */
   _completeCisTranse(product, uncleaned_fragments, frags_mapping);

   return true;
}


void ReactionEnumeratorState::_addFragments( const Molecule &fragment, 
                    const Molecule &submolecule, const Array<int> &sub_rg_atoms,
                    int *core_sub, int *core_super )
{
   const Array<int> &reactant_aam_array = _reaction.getAAMArray(_reactant_idx);
   QS_DEF(Array<int>, frag_mapping);
   
   if (!_is_frag_search)
   {
      frag_mapping.clear_resize(fragment.vertexEnd());
      frag_mapping.fffill();
      _fragments.mergeWithMolecule(fragment, &frag_mapping);
      
      /* Fragments AAM array updating */
      for (int i = 0; i < frag_mapping.size(); i++)
         _fragments_aam_array.push(0);
      for (int i = submolecule.vertexBegin(); i != submolecule.vertexEnd(); i = submolecule.vertexNext(i))
      {
         if (i < reactant_aam_array.size())
            _fragments_aam_array[frag_mapping[core_sub[i]]] = reactant_aam_array[i];
         else
            _fragments_aam_array[frag_mapping[core_sub[i]]] = 0;
      }
   }
   else
   {
      frag_mapping.clear_resize(_fragments.vertexEnd());
      for (int i = 0; i < frag_mapping.size(); i++)
         frag_mapping[i] = i;
   }
   
   QS_DEF(Array<byte>, unfrag_mon_atoms);
   unfrag_mon_atoms.clear_resize(fragment.vertexEnd());
   unfrag_mon_atoms.zerofill();
   
   _findFragAtoms(unfrag_mon_atoms, submolecule, core_sub, core_super);

   /* _is_needless array expanding */
   for (int i = _is_needless_atom.size(); i < _fragments.vertexEnd(); i++)
      _is_needless_atom.push(0);

   /* _is_needless array upadting */
   for (int i = 0; i < frag_mapping.size(); i++)
      if (unfrag_mon_atoms[i])
         _is_needless_atom[frag_mapping[i]] = 1;
   
   for (int i = 0; i < sub_rg_atoms.size(); i++)
   {
      const Vertex &sub_v = submolecule.getVertex(sub_rg_atoms[i]);
      int sub_rgroup = submolecule.getRGroups()->getSiteRGroups(sub_rg_atoms[i])[0];

      if (_att_points.size() <= sub_rgroup)
         _att_points.expand(sub_rgroup + 1);

      if (sub_v.degree() == 1)
      {
         int sub_nv_idx = sub_v.neiVertex(sub_v.neiBegin());
         int frag_rg_idx = frag_mapping[core_sub[sub_rg_atoms[i]]];
         _att_points[sub_rgroup].push(frag_rg_idx);
         if (_att_points[sub_rgroup].size() == 2)
         {
            int another_sub_v_idx = core_super[frag_mapping.find(_att_points[sub_rgroup][0])];
            if ((reactant_aam_array.size() <= another_sub_v_idx) ||
                (reactant_aam_array.size() <= sub_nv_idx))
               throw Error("Incorrect AAM");

            /* RGroup atom that have neibor with less AAM is first */
            const Vertex &another_sub_v = submolecule.getVertex(another_sub_v_idx);
            int another_sub_nv_idx = another_sub_v.neiVertex(another_sub_v.neiBegin());
            if (reactant_aam_array[another_sub_nv_idx] > reactant_aam_array[sub_nv_idx])
            {
               int tmp = _att_points[sub_rgroup][0];
               _att_points[sub_rgroup][0] = _att_points[sub_rgroup][1];
               _att_points[sub_rgroup][1] = tmp;
            }      
         }
      }
      else if (sub_v.degree() == 2)
      {
         /* ex:(N-R1-O) for R1 rgroup atom we make 2 attachment points */
         int frag_rg_idx = frag_mapping[core_sub[sub_rg_atoms[i]]];
         _att_points[sub_rgroup].push(frag_rg_idx);
         _att_points[sub_rgroup].push(frag_rg_idx);
      }
      else
         throw Error("RGroup atom can't have more than two neighbors");
   }
}

int ReactionEnumeratorState::_embeddingCallback( const Graph &subgraph, const Graph &supergraph,
                                              int *core_sub, int *core_super, void *userdata )
{
   ReactionEnumeratorState *rpe_state = (ReactionEnumeratorState *)userdata;
   const Molecule &cur_monomer = (const Molecule &)supergraph;
   const Molecule &cur_reactant = (const Molecule &)subgraph;
   QS_DEF(Molecule, submolecule);
   submolecule.clear();
   submolecule.cloneMolecule(cur_reactant, NULL, NULL);
   QS_DEF(Molecule, supermolecule);
   supermolecule.clear();
   supermolecule.cloneMolecule(cur_monomer, NULL, NULL);
   QS_DEF(Array<int>, sub_qa_array);
   sub_qa_array.clear();
   QS_DEF(Molecule, mol_fragments);
   mol_fragments.clear();

   rpe_state->_is_used_hydrogens.zerofill();

   if (!MoleculeStereocenters::checkSub(submolecule.getStereocenters(),
      supermolecule.getStereocenters(), core_sub, false))
      return 1;

   if (!MoleculeCisTrans::checkSub(submolecule, supermolecule, core_sub))
      return 1;

   Molecule &new_product = rpe_state->_full_product;


   /* Cis-Trans structure updating */
   const Array<int> &reactant_aam_array = rpe_state->_reaction.getAAMArray(rpe_state->_reactant_idx);

   QS_DEF(Array<int>, rp_mapping);
   rp_mapping.clear_resize(submolecule.vertexEnd());
   rp_mapping.fffill();

   rpe_state->_findR2PMapping(submolecule, new_product, rp_mapping);

   rpe_state->_cistransUpdate(submolecule, supermolecule, NULL, rp_mapping, core_sub);
   rpe_state->_stereocentersUpdate(submolecule, supermolecule, rp_mapping, core_sub, core_super);


   /* Finding indices of ignored query atoms */
   for (int i = submolecule.vertexBegin(); i != submolecule.vertexEnd(); i = submolecule.vertexNext(i))
   {
      if (submolecule.getQueryAtom(i).type == QUERY_ATOM_RGROUP)
         sub_qa_array.push(i);
   }

   mol_fragments.cloneMolecule(supermolecule, NULL, NULL);

   /* Updating fragments molecule */
   rpe_state->_addFragments(mol_fragments, submolecule, sub_qa_array, core_sub, core_super);
   
   MoleculeRGroups *sub_rgroups = submolecule.getRGroups();

   if (rpe_state->is_self_react)
   {
      ReactionEnumeratorState self_rxn_rpe_state(*rpe_state);
      self_rxn_rpe_state._is_frag_search = true;
      self_rxn_rpe_state._reactant_idx = self_rxn_rpe_state._reaction.reactantNext(self_rxn_rpe_state._reactant_idx);
      
      if (self_rxn_rpe_state._reactant_idx != self_rxn_rpe_state._reaction.reactantEnd())
         self_rxn_rpe_state._start_ee(mol_fragments);
   }

   ReactionEnumeratorState new_rpe_state(*rpe_state);
   
   new_rpe_state._reactant_idx = new_rpe_state._reaction.reactantNext(new_rpe_state._reactant_idx);
   new_rpe_state.buildProduct();
   return 0;
}